rm(list = ls())

#setwd("YOUR-WORKING-DIRECTORY")

# install.packages(c("doParallel", "foreach", "magic"))

library("doParallel")
library("foreach")
library("magic")
library("MASS")
library("beepr")
library("list")

set.seed(88)

# stopCluster(cl)
# cl <- makeCluster(detectCores() - 2)
# registerDoParallel(cl)

source("replication-functions.R")


finished <- read.csv("ictreg-fig-1.csv")

# overall.1000 <- subset(finished, case == "overall" & n == 1000)
# overall.2500 <- subset(finished, case == "overall" & n == 2500)
# overall.1e04 <- subset(finished, case == "overall" & n == 1e04)

subgroup.1000 <- subset(finished, case == "subgroup" & n == 1000)
subgroup.2500 <- subset(finished, case == "subgroup" & n == 2500)
subgroup.1e04 <- subset(finished, case == "subgroup" & n == 1e04)

corr1.1000 <- subset(finished, case == "corr1" & n == 1000)
corr1.2500 <- subset(finished, case == "corr1" & n == 2500)
corr1.1e04 <- subset(finished, case == "corr1" & n == 1e04)

corr2.1000 <- subset(finished, case == "corr2" & n == 1000)
corr2.2500 <- subset(finished, case == "corr2" & n == 2500)
corr2.1e04 <- subset(finished, case == "corr2" & n == 1e04)

corr3.1000 <- subset(finished, case == "corr3" & n == 1000)
corr3.2500 <- subset(finished, case == "corr3" & n == 2500)
corr3.1e04 <- subset(finished, case == "corr3" & n == 1e04)

n.bootstrap <- 5000

# corrs <- c(1.0, 0.0, 0.0, 
# 	       0.0, 1.0, 0.0, 
# 	       0.0, 0.0, 1.0)


# pop.data <- generateComparisonData(n = 1e07, corrs = corrs, gammas = rep(0.5, 3), 
# 	deltas = rep(0.5, 3), n.control.items = 4, p.t = 0.5, n.groups = 5)


# {
# 	par(cex = 0.7)
# 	hist(pop.data$x1[pop.data$group == 1], xlim = c(-7, 7), col = "black", freq = FALSE, 
# 		ylim = c(0, 0.5))
# 	shades <- seq(1, 0, length.out = 4)
# 	for (i in 2:5) {
# 		hist(pop.data$x1[pop.data$group == i], 
# 			col = rgb(0, 0, 1, shades[i - 1]), add = TRUE, freq = FALSE)
# 	}
# }

# necFunc <- c("adiag", "coef.ictreg", "vcov.ictreg", "ginv")


# # # case 1
# # {

# # 	h <- mean(pop.data$z)

# # 	overall.1000 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# # 	  generateComparison(n.sample = 1000, data = pop.data, n.control.items = 4, 
# # 		h = h, grouping = "overall")	

# # 	}


# # 	overall.2500 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# # 	  generateComparison(n.sample = 2500, data = pop.data, n.control.items = 4, 
# # 		h = h, grouping = "overall")	

# # 	}


# # 	overall.1e04 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# # 	  generateComparison(n.sample = 1e04, data = pop.data, n.control.items = 4, 
# # 		h = h, grouping = "overall")	

# # 	}

# # }


# # case 2
# { 

# 	h <- tapply(pop.data$z, pop.data$group, mean)

# 	subgroup.1000 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparison(n.sample = 1000, data = pop.data, n.control.items = 4, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}


# 	subgroup.2500 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparison(n.sample = 2500, data = pop.data, n.control.items = 4, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}


# 	subgroup.1e04 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparison(n.sample = 1e04, data = pop.data, n.control.items = 4, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}

# }


# # case 3
# {

# 	corrs <- c(1.0, 0.0, 0.5, 
# 		       0.0, 1.0, 0.0, 
# 		       0.5, 0.0, 1.0)

# 	pop.data <- generateComparisonData(n = 1e07, corrs = corrs, gammas = rep(0.5, 3), 
# 		deltas = rep(0.5, 3), n.control.items = 4, p.t = 0.5, n.groups = 5)

# 	h <- tapply(pop.data$z, pop.data$group, mean)

# 	corr1.1000 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparison(n.sample = 1000, data = pop.data, n.control.items = 4, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}

# 	corr1.2500 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparison(n.sample = 2500, data = pop.data, n.control.items = 4, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}


# 	corr1.1e04 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparison(n.sample = 1e04, data = pop.data, n.control.items = 4, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}

# }


# # case 4
# { 

# 	corrs <- c(1.0, 0.0, 0.5, 
# 		       0.0, 1.0, 0.5, 
# 		       0.5, 0.5, 1.0)

# 	pop.data <- generateComparisonData(n = 1e07, corrs = corrs, gammas = rep(0.5, 3), 
# 		deltas = rep(0.5, 3), n.control.items = 4, p.t = 0.5, n.groups = 5)

# 	h <- tapply(pop.data$z, pop.data$group, mean)

# 	corr2.1000 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparison(n.sample = 1000, data = pop.data, n.control.items = 4, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}

# 	corr2.2500 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparison(n.sample = 2500, data = pop.data, n.control.items = 4, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}


# 	corr2.1e04 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparison(n.sample = 1e04, data = pop.data, n.control.items = 4, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}

# }


# # case 5
# { 

# 	corrs <- c(1.0, 0.5, 0.5, 
# 		       0.5, 1.0, 0.5, 
# 		       0.5, 0.5, 1.0)

# 	pop.data <- generateComparisonData(n = 1e07, corrs = corrs, gammas = rep(0.5, 3), 
# 		deltas = rep(0.5, 3), n.control.items = 4, p.t = 0.5, n.groups = 5)

# 	h <- tapply(pop.data$z, pop.data$group, mean)

# 	corr3.1000 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparison(n.sample = 1000, data = pop.data, n.control.items = 4, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}

# 	corr3.2500 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparison(n.sample = 2500, data = pop.data, n.control.items = 4, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}


# 	corr3.1e04 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparison(n.sample = 1e04, data = pop.data, n.control.items = 4, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}

# }

# beep()

pdf(file = "figure-4.pdf", width = 10, height = 10)

{


index <- rep(0:1, n.bootstrap)

par(mfrow = c(3, 3), cex = 0.75, oma = c(2, 2, 2, 12), mar = c(5, 3, 2, 2) + 0.05, xpd = NA)
	# bias
	for (i in 1:3) {

		plot(x = 1:3, y = c(mean(subgroup.1000[index == 0, i]), 
			mean(subgroup.2500[index == 0, i]), 
			mean(subgroup.1e04[index == 0, i])
			), 
	        lwd = 1.5, 
		    type = "b", 
		    ylim = c(0.485, 0.525), 
		    yaxt = "n",
		    ylab = "Bias", 
		    xaxt = "n", 
		    xlab = "Sample Size", 
		    main = substitute(delta[m], list(m = i - 1)), 
		    lty = 1, 
		    pch = 1
		)

		par(xpd = FALSE); abline(h = 0.5, col = "darkgray"); par(xpd = NA)
	    axis(side = 2, at = seq(0.49, 0.52, 0.01), labels = seq(0.49, 0.52, 0.01) - .5)
	    axis(side = 1, at = 1:3, labels = c("1000", "2500", "10000"))

		lines(x = 1:3, y = c(mean(subgroup.1000[index == 1, i]), 
			mean(subgroup.2500[index == 1, i]), 
			mean(subgroup.1e04[index == 1, i])
			), 
		    lwd = 1.5, 
		    type = "b", 
		    col = "black", 
		    lty = 1, 
		    pch = 16
		   )

		lines(x = 1:3, y = c(mean(corr1.1000[index == 1, i]), 
			mean(corr1.2500[index == 1, i]), 
			mean(corr1.1e04[index == 1, i])
			), 
	        lwd = 1.5, 
		    type = "b", 
		    col = "red", 
		    lty = 2, 
		    pch = 2
		   )

		lines(x = 1:3, y = c(mean(corr2.1000[index == 1, i]), 
			mean(corr2.2500[index == 1, i]), 
			mean(corr2.1e04[index == 1, i])
			), 
	        lwd = 1.5, 
		    type = "b", 
		    col = "blue", 
		    lty = 2, 
		    pch = 17
		   )

		lines(x = 1:3, y = c(mean(corr3.1000[index == 1, i]), 
			mean(corr3.2500[index == 1, i]), 
			mean(corr3.1e04[index == 1, i])
			), 
	        lwd = 1.5, 
		    type = "b", 
		    col = "darkgreen", 
		    lty = 2, 
		    pch = 4
		   )

	}


	legend("topright", legend = c("No aux. info",  
	  "No segregation", "Segregation on X1", "Segregation on X1, X2", "X1, X2 correlated"), 
	  col = c("black", "black", "red", "blue", "darkgreen"), 
	  lty = c(1, 1, 1, 2, 2),
	  pch = c(1, 16, 2, 17, 4),  
	  lwd = 1.5, 
	  bty = "n", 
	  seg.len = 1.5, 
	  inset = c(-1, 0)
	  )


	# rmse
	for (i in 1:3) {

		plot(x = 1:3, y = c(rmse(subgroup.1000[index == 0, i], 0.5), 
			rmse(subgroup.2500[index == 0, i], 0.5), 
			rmse(subgroup.1e04[index == 0, i], 0.5)
			), 
		    lwd = 1.5, 
		    type = "b", 
		    ylim = c(0, 0.2), 
		    ylab = "RMSE", 
		    xaxt = "n", 
		    xlab = "Sample Size", 
		    main = substitute(delta[m], list(m = i - 1)), 
		    lty = 1, 
		    pch = 1
		)
	    axis(side = 1, at = 1:3, labels = c("1000", "2500", "10000"))

		par(xpd = FALSE); abline(h = 0, col = "darkgray"); par(xpd = NA)

		lines(x = 1:3, y = c(rmse(subgroup.1000[index == 1, i], 0.5), 
			rmse(subgroup.2500[index == 1, i], 0.5), 
			rmse(subgroup.1e04[index == 1, i], 0.5)
			), 
		    type = "b", 
	        lwd = 1.5, 
		    col = "black", 
		    lty = 1, 
		    pch = 16
		   )

		lines(x = 1:3, y = c(rmse(corr1.1000[index == 1, i], 0.5), 
			rmse(corr1.2500[index == 1, i], 0.5), 
			rmse(corr1.1e04[index == 1, i], 0.5)
			), 
		    type = "b", 
	        lwd = 1.5, 
		    col = "red", 
		    lty = 2, 
		    pch = 2
		   )

		lines(x = 1:3, y = c(rmse(corr2.1000[index == 1, i], 0.5), 
			rmse(corr2.2500[index == 1, i], 0.5), 
			rmse(corr2.1e04[index == 1, i], 0.5)
			), 
		    type = "b", 
	        lwd = 1.5, 
		    col = "blue", 
		    lty = 2, 
		    pch = 17
		   )

		lines(x = 1:3, y = c(rmse(corr3.1000[index == 1, i], 0.5), 
			rmse(corr3.2500[index == 1, i], 0.5), 
			rmse(corr3.1e04[index == 1, i], 0.5)
			), 
		    type = "b", 
	        lwd = 1.5, 
		    col = "darkgreen", 
		    lty = 2, 
		    pch = 4
		   )

	}

	legend("topright", legend = c("No aux. info",  
	  "No segregation", "Segregation on X1", "Segregation on X1, X2", "X1, X2 correlated"), 
	  col = c("black", "black", "red", "blue", "darkgreen"), 
	  lty = c(1, 1, 1, 2, 2),
	  pch = c(1, 16, 2, 17, 4), 
	  lwd = 1.5, 
	  bty = "n", 
	  seg.len = 1.5, 
	  inset = c(-1, 0)
	  )

	# coverage
	for (i in 1:3) {

		plot(x = 1:3, y = c(
				cover(x   = subgroup.1000[index == 0, i], 
					  se  = subgroup.1000[index == 0, i + 6], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = subgroup.2500[index == 0, i], 
					  se  = subgroup.2500[index == 0, i + 6], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = subgroup.1e04[index == 0, i], 
					  se  = subgroup.1e04[index == 0, i + 6], 
					  par = 0.5, 
					  alpha = 0.05)
				), 
		    lwd = 1.5, 
		    type = "b", 
		    ylim = c(0.9, 1), 
		    ylab = "Coverage", 
		    xaxt = "n", 
		    xlab = "Sample Size", 
		    main = substitute(delta[m], list(m = i - 1)), 
		    lty = 1, 
		    pch = 1
		)
		
	    axis(side = 1, at = 1:3, labels = c("1000", "2500", "10000"))
		par(xpd = FALSE); abline(h = 0.95, col = "darkgray"); par(xpd = NA)

		lines(x = 1:3, y = c(
				cover(x   = subgroup.1000[index == 1, i], 
					  se  = subgroup.1000[index == 1, i + 6], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = subgroup.2500[index == 1, i], 
					  se  = subgroup.2500[index == 1, i + 6], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = subgroup.1e04[index == 1, i], 
					  se  = subgroup.1e04[index == 1, i + 6], 
					  par = 0.5, 
					  alpha = 0.05)
				), 
	        lwd = 1.5, 
		    type = "b", 
		    lty = 1, 
		    pch = 16
		   )

		lines(x = 1:3, y = c(
				cover(x   = corr1.1000[index == 1, i], 
					  se  = corr1.1000[index == 1, i + 6], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = corr1.2500[index == 1, i], 
					  se  = corr1.2500[index == 1, i + 6], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = corr1.1e04[index == 1, i], 
					  se  = corr1.1e04[index == 1, i + 6], 
					  par = 0.5, 
					  alpha = 0.05)
				), 
	        lwd = 1.5, 
		    type = "b", 
		    col = "red", 
		    lty = 2, 
		    pch = 2
		   )

		lines(x = 1:3, y = c(
				cover(x   = corr2.1000[index == 1, i], 
					  se  = corr2.1000[index == 1, i + 6], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = corr2.2500[index == 1, i], 
					  se  = corr2.2500[index == 1, i + 6], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = corr2.1e04[index == 1, i], 
					  se  = corr2.1e04[index == 1, i + 6], 
					  par = 0.5, 
					  alpha = 0.05)
				), 
	        lwd = 1.5, 
		    type = "b", 
		    col = "blue", 
		    lty = 2, 
		    pch = 17
		   )


		lines(x = 1:3, y = c(
				cover(x   = corr3.1000[index == 1, i], 
					  se  = corr3.1000[index == 1, i + 6], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = corr3.2500[index == 1, i], 
					  se  = corr3.2500[index == 1, i + 6], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = corr3.1e04[index == 1, i], 
					  se  = corr3.1e04[index == 1, i + 6], 
					  par = 0.5, 
					  alpha = 0.05)
				), 
	        lwd = 1.5, 
		    type = "b", 
		    col = "darkgreen", 
		    lty = 2, 
		    pch = 4
		   )

		}

		legend("topright", legend = c("No aux. info",  
		  "No segregation", "Segregation on X1", "Segregation on X1, X2", "X1, X2 correlated"), 
		  col = c("black", "black", "red", "blue", "darkgreen"), 
		  lty = c(1, 1, 2, 2, 2),
		  pch = c(1, 16, 2, 17, 4), 
		  lwd = 1.5, 
		  bty = "n", 
		  seg.len = 1.5, 
		  inset = c(-1, 0)
		  )

  # title(main = c("Statistical Properties of the Augmented List Experiment Estimator"), outer = TRUE)


    

}

dev.off()




# finished <- rbind(
# 	# cbind(overall.1000, n = 1000, case = "overall"), 
# 	# cbind(overall.2500, n = 2500, case = "overall"), 
#     # cbind(overall.1e04, n = 1e04, case = "overall"), 

# 	cbind(subgroup.1000, n = 1000, case = "subgroup"), 
# 	cbind(subgroup.2500, n = 2500, case = "subgroup"), 
# 	cbind(subgroup.1e04, n = 1e04, case = "subgroup"), 

# 	cbind(corr1.1000, n = 1000, case = "corr1"), 
# 	cbind(corr1.2500, n = 2500, case = "corr1"), 
# 	cbind(corr1.1e04, n = 1e04, case = "corr1"), 

# 	cbind(corr2.1000, n = 1000, case = "corr2"), 
# 	cbind(corr2.2500, n = 2500, case = "corr2"), 
# 	cbind(corr2.1e04, n = 1e04, case = "corr2"), 

# 	cbind(corr3.1000, n = 1000, case = "corr3"), 
# 	cbind(corr3.2500, n = 2500, case = "corr3"), 
# 	cbind(corr3.1e04, n = 1e04, case = "corr3")
# 	)

# finished <- as.data.frame(finished)
# names(finished)[13] <- "augmented"

# write.csv(finished, "ictreg-fig-1.csv", row.names = FALSE)


# stopCluster(cl)




# # q("no")